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Abstract 

We investigate a hybrid PDE/Monte Carlo technique for the variance reduced simulation 
of an agent-based multiscale model for tumor growth. The variance reduction is achieved by 
combining a simulation of the stochastic agent-based model on the microscopic scale with 
a deterministic solution of a simplified (coarse) partial differential equation (PDE) on the 
macroscopic scale as a control variable. We show that this technique is able to significantly 
reduce the variance with only the (limited) additional computational cost associated with 
the deterministic solution of the coarse PDE. We illustrate the performance with numerical 
experiments in different regimes, both in the avascular and vascular stage of tumor growth. 


1 Introduction 


Tumor growth is a complex biological phenomenon consisting of processes on different scales. 
On the cellular level - which will be referred to as the microscopic scale in this paper - one 
has to track the random motion of cells, as well as the cell division and cell death. The latter 
are governed by numerous intracellular processes. Furthermore, the cellular behavior is strongly 
coupled to the environment and vice versa. For example, cell proliferation is determined by 
the local oxygen concentration and the local cell density while hypoxia on the other hand can 
trigger apoptosis, but cells also consume oxygen. This two-way feedback creates a very specific 
dynamics characterizing the development of the tumor. A hypoxic zone develops in the middle 
of the tumor, which in turn triggers endothelial cells to vascularize the tumor. This process, 
also known as angiogenesis |^, ensures that the tumor’s need for oxygen and other nutrients is 
satisfied, which implies that the tumor can grow further. 

Smaller avascular tumors can be easily simulated on the microscopic scale using agent-based 
models. We can distinguish two classes of models. On one hand, cellular automata update grid 
cells based on a number of phenomenological rules 23 27 , while on the other hand lattice-free 


models typically consist of a set of ordinary differential equations (ODEs) attached to each cell. 
On long time scales, we are typically interested in the tumor as a whole, which we call the macro¬ 
scopic scale. Agent-based models are typically not well suited to use on this larger scales, since 
the individual based character implies a large computational cost for a large number of particles, 
corresponding to larger tumors. One may choose to model the system directly on this scale using 
continuum models, based on mass balance equations [3l[^[M|[^[Ml[^ . While this approach is 


*Department of Computer Science, KU Leuven, Celestijnenlaan 200A, 3001 Leuven, Belgium (first- 
name.lastname@cs.kuleuven.be). The first author’s work was supported by the Agency for Innovation by Science 
and Technology in Flanders (IWT) 


1 






significantly cheaper and easier to analyze than agent-based models, it cannot capture discrete 
features as branching of a vascular network or events regulated by intracellular concentrations. 
This insight gave rise to multiscale models where agent-based models are typically used to model 
the cellular component, while the environment is mostly described by a set of reaction-diffusion 
partial differential equations (PDEs), corresponding to the macroscopic scale. Examples can be 
found in fTl[8 17 28 29 . For a review about the current state of the art in multiscale-modeling 


of tumor growth, we refer to |^. 


Due to the random motion and the influence on the environment, the simulations are subject 
to noise. When simulating with a standard Monte Carlo algorithm, the variance can only be 
reduced by increasing the number of particles at the cost of computational efficiently. Various 
techniques for variance reduction such as antithetic variables, control variates and importance 
sampling are described in literature, see e.g. |4 24 for an overview. Recently, several hybrid 
PDE/Monte Carlo algorithms have been proposed in the literature to achieve variance reduction 
by coupling a PDE-based discretization to a Monte Carlo simulation [To|[^[^ . 

The contributions of this paper are two-fold: 

• We develop a multiscale model where the random motion is modeled using stochastic 
differential equations (SDEs), the intracellular variables for the cell cycle and apoptosis are 
described by ODEs and the environment, consisting of diffusible components, is modeled 
by PDEs. The model is a modified version of the cellular automaton model of Owen et 
al. 29 . The main differences are that the new model is lattice-free and the fact that our 


model does not contain any explicit delay terms. 


• We propose a novel technique to reduce the variance on the results of the adapted multi¬ 
scale agent-based model. Based on the ideas in [^[^, we develop a hybrid PDE/Monte 
Carlo method using a coarse stochastic process (called the control process) and a corre¬ 
sponding PDE. In this specihc case, the control process modeling the spatial behavior of 
the individual cells contains all details of the microscopic model except for cell births, cell 
deaths and VEGF secretion. The keypoint is to obtain this missing information with re¬ 
duced variance by an appropriate coupling between the full microscopic agent-based model 
and the control process. 


We first give a detailed overview of the different layers of the model. Next, we describe the 
variance reduction algorithm in the section]^ We illustrate the technique numerically in section]^ 
Finally, in section we elaborate on a few possibilities for future research. 


2 Models 


In this section, we describe a multiscale model for tumor growth. The microscopic model for 
tumor growth is based on the ideas used to describe bacterial chemotaxis 11,32 , the multiscale 


cellular automaton model was developed by Owen and coworkers 29 and our goal is to perform 


variance reduction in order to estimate the resulting population densities in a more accurate way. 
As in 32 , the model proposed in this paper is time and space continuous. Apart from the fact 


that the model is lattice-free, making the computational cost quasi independent of the size of 
the domain and hence, we can easily rescale the system to simulate larger tumors (compared to 
the examples given in [^). 

We distinguish two main components: the environment, modeled by a couple of reaction diffu¬ 
sion equations and the agent-based model describing the individual cellular motion and internal 
variables (e.g. cell cycle, apoptosis state and internal concentrations such as VEGF and p53) 
attached to each cell. 
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We consider three types of cells, indexed by 1 < p < P = 3: normal cells {p = 1), cancer cells 
(p = 2), and endothelial cells (that build up blood vessels, p = 3). For each of these cell types, 
we consider an ensemble of Ip{t) cells, and consider three state variables: position a: S cell 
cycle phase cj) G [0,1]. The intracellular concentrations [VEGFjint, [p53] and apoptosis variable z 
are scalars S M. These cells evolve according to evolution laws that depend on the concentration 
[02\{x,t) of oxygen and [VEGF](a;,<) of the Vascular Endothelial Growth Eactor (which we call 
the environment). 

Remark that the reaction-diffusion partial differential equations (PDEs) describing the diffusible 
components of the environment still need to be solved on a grid, but this cost is marginal due to 
the sparsity of the involved linear systems, which ensures that the cost dependent on the domain 
size is limited. 

We now give an overview of the notations that will be used throughout the paper, after which 
we describe the evolution laws for the environment, and detail the evolution laws for each of the 
cell types. The cell type dependency is mainly caused by cell type dependent coefficients, which 
will be discussed later on in the section describing the agent-based model in more detail. 

Notation 

• The state variables attached to a single cell of type p at time t are position Xp{t), cell 

cycle phase generation Cpit): internal concentrations [p53]p(t), [VEGEint]p(t) and 

apoptosis variable Zp{t). 

• Particle number densities are denoted by np{x, t) indexed by a suitable subscript to indicate 
the nature of the density. Eurther, ny{x,t) is used to describe the vascular density. 

• To keep a consistent notation throughout the paper, we introduce the following convention. 
If, at a moment t = G, the cell with index i* in population p divides, we set 

ipin = 4(G) +1 (1) 

in which the symbol t*_ is used to emphasize that the involved number of cells is meant 
to be taken just before the division. Simultaneously, we introduce a new cell as specified 
in the paragraph concerning cell division (see page|^. When a cell undergoes apoptosis, 
it is removed from the simulation. To avoid cumbersome renumbering of the cells in the 
text, we associate a weight Wi^p{t) to each of the cells. If the cell is alive, the corresponding 
weight is one; upon apoptosis, it becomes zero. The active number of cells is therefore: 

Ip (t) 

Ipit) = y ' Wi,p(t) ( 2 ) 

i=l 

• The evolution of the state variables is influenced in various ways by the (local) environment. 
The latter will be modelled by means of diffusible components [VEGE](a:,t) describing the 
VEGE concentration, while [02]{x,t) denotes the oxygen concentration. 

2.1 Agent-based model 

In this section we give a detailed overview of the evolution of the different state variables attached 
to each cell of the different cellular populations. 
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Position. The random motion of the position of the cells is described as a biased Brownian 
motion. Cells of type p move randomly with diffusion coefficient Dp, and the cells are chemotac- 
tically attracted towards high concentrations with sensitivity Xp- This sensitivity is especially 
important for the endothelial cells, responsible for blood vessel growth, 

dXp{t) = Xpy[VEGF]{Xp{t),t)) (l- 2 Ei^pMhR\ at + v^dfTt (3) 

\ ^max,p / 

The cell number density can be computed as: 

/p(i) 

np{x,t) = (4) 

i=l 


where Ip{t) denotes the total number of cells of type p at time t and S denotes the classical 
Dirac kernel, resulting in a standard histogram. Remark, in contrast to the cellular automaton 
model described in 29 , the resulting equation for the position is a stochastic differential equation 
(SDE) instead of the discrete space jumps used in 29 . Finally, we have to stress the fact that 


the above equation is general for all the cell types. To be more concrete, normal cells don’t 
move at all, while cancer cells are characterized by pure diffusive motion and endothelial cells 
demonstrate diffusive behavior but they also respond to chemotactic cues. 


Cell division is modeled by means of the following ODE: 


d$p(t) 

dt 


[02]iXp{t),t) 


+ [02\iXp(t), t)) 


H(Cp(t)-Cp.r 


(5) 


where 


min,p 


denotes the minimal time needed for a cell to complete one cell cycle and ^ 


indicates the generation of a cell. Remark that Tmin.p depends on the cell type. To be more 
specific, cancer cells are able to proceed twice as fast as normal cells during the cell cycle in 
a given environment (see table [^. Naturally the cell cycle speed depends on the local oxygen 
concentration [02]{Xp{t),t) as observed by the cell while evolving through the cycle. The higher 
the oxygen concentration, the faster the cycle is completed, while the cell cycle is put on hold 
when the cell suffers from hypoxia. A more detailed biological motivation for this model can be 
found in and its supplementary material. Remark that in the cellular automaton model 

by Owen et al (see [^) all cells can divide an unlimited number of times, which corresponds to 
the hypothesis that all cells are stem cells, which is obviously not a realistic assumption. Thus, 
we have extended the model to account for the fact that cells are only able to divide a finite 
number of times (i.e. Cp,max)- To be more specific we added a factor H(^p(t) — Cp.max) to check 
for the generation of the corresponding cells. Here, we assume that normal cells can divide only 


4 times, which is consistent with 12 . On the other hand we assume that all the cancer cells are 


cancer stem cells, which is still a simplification. 

If, for the cell with index i* in population p at time t = t*, we obtain <i)(t*) > 1, we introduce a 
new cell in the simulation. We adjust Ip{t) according to equation 0 and set = 0. and 

the generation of the parent cell increases by one. The new cell inherits the state from the cell 
that divides except for the generation (: 


Xj^aipA) = AxpA) 
[p53]p^p),p(t*) = [p53],*,p(t*) 


Eip{t),p{t ) — Zi*^p{t ) 

[VEGFi„t]p^(t),p(t*) = [VEGFi„t].*.p(t*) 

C/,(t).p(t*) = 0 


( 6 ) 
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Intracellular model. We introduce a intracellular module consistent with 29 in order to 


describe some important intracellular concentrations, namely the p53 concentration [p53] and 
the intracellular VEGF concentration [VEGFint]- The former can be seen as an estimator for 
the number of mutations that a cell has undergone during its lifetime. We have: 


d[p53]p(t) 

dt 


= Ci-C2 


Cp53 + [02]{Xp{t),t) 


[p53]p(t) 


d[VEGFi„t]p(t) [p53]p(t)[VEGFi„t]p(t) 

dt """ Js + [VEGFi„t]p(t) 


[02]{Xp{t),t) 

^CveGF + [02]{Xp{t),t) 


[VEGFi„t]p(t) 


(7) 


Cells are storing VEGF intracellular (i.e. [VEGFjint) during hypoxic conditions and release it 
once this intracellular concentration has reached a certain threshold level [VEGFi„t]thr- Further, 
Cl, ... C5 and Cpsa, CVegf are constants that can be found in table[^ Next we describe the model 
for apoptosis, depending on the cell type. Therefore we formally define 7apt.p(zi ^^-p) = Fp{z,np) 
as the apoptosis rate, which is further specified in the following paragraphs. 


Apoptosis for normal cells. For normal cells, cell death is completely determined by the 
subcellular p53-concentration. So, we set the apoptosis variable z := [p53]. The apoptosis 
threshold 7apt can then then be written as: 

ls.pt, i{z,ni) = H (z - ZhighH(nthr - «i) - ziowH(ni - nthr)) (8) 

where H indicates the Heaviside function. This definition of 7apt implies that normal cells undergo 
apoptosis if 7apt('2, ni) = 1 corresponding to the situation that z has reached a certain threshold 
value depending on the harshness of the environment. The threshold value is lower in case of a 
harsh environment, defined as ni < nthr, where nthr denotes a threshold value for the normal 
cells. 


Apoptosis for cancer cells. In contrast to normal cells, the apoptosis mechanism for tumor 
cells is independent of the p53-concentration since this mechanism to regulate the normal cell 
cycle does not function properly anymore in a tumor. Gancer cells are able to go into a quiescent 
state when expressed to hypoxic circumstances, meaning that they don’t consume any nutrients 
anymore for a while. However the duration of this quiescent state is limited, which implies that 
cancer cells will also undergo apoptosis when the hypoxia holds too long. On the other hand, 
cancer cells have the ability to recover quickly once there is again more oxygen available. This 
mechanism can be modeled by the following equation: 

= AH([ 02 ]thr - [02]{Xp{t),t))- BZ{t)li{[02]{Xp{t),t) - [ 02 ]thr) 

at - V -^ ^- V -^ 

Linear increase during hypoxia Exponential decay if [02](Xp(f),t)>[02lthr 

where A, B are constants. Further, the first term models the hypoxic state, i.e. the local oxygen 
concentration [02]{Xp{t),t) drops below the threshold level [02]thr- During this hypoxic period, 
the internal variable z increases steadily. On the other hand, the second term describes the 
recovery of the cancer cells if the environment is not hypoxic anymore, which is captured by the 
exponential decay term of Z(t). Gancer cells die if Z(t) > 1, corresponding to 7 apt, 2 (-z) = H(z—1). 
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Endothelial Cells. Remark that the model equations concerning cell division and cell death 
will not be used for endothelial cells. Consistent with existing literature, the so-called snail-trail 
approach is used to model sprouting angiogenesis [2 19 . 


Full agent-based model. This results in the following set of equations for the full agent-based 
model: 


dXp{t) = Xp^[VEGF]{Xp{t),t)) 1 - 


d^p(t) _ 

dt 

d[p53]p(f) 
dt 

d[VEGFint]p(t) 

dt 


[02]iXpit),t) 


,p{C4>,p + [02]{Xp{t),t)) 
[02]{Xp{t),t) 


^max,p / 

H iCpi'l) ~ Cp ,max) 


= Cl - C2 


Cp53 + [02]{Xp(t),t) 


[p53]p(t) 


[p53]p(t)[VEGFi„t]p(t) 
J 5 + [VEGFi„4(t) 


[02]{Xpit),t) 


Gvegf + [02]iXp(t),t) 


[VEGFi„t]p(t) 


.qapt,p(^5 t) — Fpi^Z^Hp^t') 

The corresponding parameter values can be found in table 


(9) 


Parameter 

ni 

n2 

n3 

units 

Xp 

0.0 

0.0 

2 X 10-^ 

cm2/min/nM 

Cp,niax 

4 

CX) 

4 

times 

C^^p 

3 

1.4 


mmHg 

CvEGF 

0.01 

0.01 

0.01 

mmHg 

Cp53 

0.01 

0.01 

0.01 

mmHg 

^p,min 

3000 

1600 


min 

-^-high 

0.8 



dimensionless 

•^low 

0.08 



dimensionless 

?^thr 

0.75 



dimensionless 

[02]thr 


8.9 


mmHg 

[VEGFi„t]thr 

0.27 

0.27 


nM 

Cl 

2 X 10-3 

2 X 10-3 


min'^ 

C2 

1 X 10-2 

1 X 10-2 


min'^ 

C 3 

2 X 10-3 

2 X 10-3 


min'^ 

C4 

2 X 10-3 

2 X 10-3 


min'^ 

C 5 

1 X 10-2 

1 X 10-2 


min'^ 

J 5 

0.04 

0.04 


nM 

A 


1 


min'^ 

B 


2.5 X 10-3 


min'^ 


Table 1: Parameter values related to the populations. 


2.2 Coarse Description 

An alternative approach to model tumor growth is to describe the evolution of the populations 
as a whole in a probabilistic way using partial differential equations (PDFs). In general, this 
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approach yields a reaction-diffusion PDE. However, in this case it is not possible to derive a 
closed formulation for the reaction terms since those all depend on intracellular variables. A 
continuum description for the model outlined above without births and deaths can be found 
in 20 . The resulting macroscopic equation - achieved by taking the limit for a high number of 


particles - for the evolution of the populations reads: 


dtTip{x^t^ — Dp^ n,p{x^t^ • 




- '^[VEGF]{x, t) 

\ ^o.max / 


( 10 ) 


where no reactions (cell divisions, cell deaths) are taken into account. Next, we introduce the 
following macroscopic time-stepper: 


np{x, t^+^) = np{x, t^) -\- 5tDpS/^np{x, t^) 


- 5t XpV. 


np{x,t^) ( 1 - j \/[VEGF]{x,t’^'' 

^p,max 


( 11 ) 


which uses a first order Euler discretization to discretize the time derivative and a second order 
central finite volume scheme to discretize the spatial derivatives. Further details can be found 
in section |4l 


2.3 Angiogenesis 

The growth of new blood vessels, also known as angiogenesis is essential for the development of 
a tumor. Hanahan and co-authors identified it as one of the hallmarks of cancer (see 16 Ell) . In 
the early stages of cancer, the existing vasculature is able to provide enough oxygen and other 
nutrients. But as soon as the size of the tumor has reached a certain threshold, a hypoxic zone 
develops in the middle of the tumor. To cope with this phenomenon, the tumor secretes VEGF, 
a growth factor, which triggers endothelial cells to move chemotactically towards this hypoxic 
zone and grow new blood vessels. In this paper we will use an existing model for angiogenesis, 
described in 29 . We distinguish two phenotypes: endothelial cells can either be motile leader 


cells (also called tip cells) or static stalk cells. We model proliferation of endothelial cells by means 
of the so-called snail-trail approach, where each tip cell produces a new (static) endothelial cell at 
its previous position, creating a trail of static stalk cells behind him. Apart from this feature, new 
tip cells -known as sprouts- can emerge from active vessels with sprouting probability Psprout 
along the active vessels, (see 29 ): 


„ ,, Pr^.4VEGF]ix,t) 

^P^o-‘-^V.p,out + [I^AGP](x,t) 


( 12 ) 


where Pmax = 3 x 10 ^ min“^ indicates the maximal endothelial sprouting rate ( see 29 ) and 


I^sprout= 0.5 nM denotes the VEGF concentration at which the sprouting probability is half 
maximal. Remark that the probability of the emergence of two sprouts close to each other 
within the same timestep is zero. A biological explanation for this fact can be found in the 

and (^[^[22 


supplementary material provided with 29 


where the authors pointed out that 
delta-notch signaling inhibits the formation of new sprouts in neighbouring endothelial cells. 
Additionally, the vessel radii are also adapted dynamically based on the work of 29 where 


pruning of the vessels was also incorporated in the model: if the pressure in a branch is too low, 
the corresponding will collapse. 
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2.4 Environment 


The cellular environment consists of two diffusible components regulating the behavior of the 
cells in various ways. Oxygen is evidently important for the cells to proceed through the cell 
cycle. The local oxygen concentration is determined from the following equation: 

dt[02]ix,t) = D[02]'^^[02]{x,t) + il^[02]n.u{x,t){[02]h\ood - [02]ix,t)) 

'• -V-^ '-V-' 

diffusion exchange with blood 

(13) 

+ [02]{x,t]k[02]^np{x,t) 

p^l 

^ "V 

Consumption 


where T^[ 02 ] Hi® diffusion coefficient of oxygen, 'ipl 02 ] denotes the permeability of the oxygen 
through the vessels. ny{x,t) describes the surface area occupied by the vessel at position x. 
[O 2 ]biood(a;,i) = [ 02 ]reiH{x,t)/Hin defines the oxygen concentration in a blood vessel located at 
position X. [02]ref is a reference oxygen concentration, H{x,t) is the haematocrit at location x 
and time t and Hi^ is the haematocrit at an inflow node and by default it is set to Hi^ — 0 . 45 . 
The last term in reflects the fact that all cells consume oxygen at a rate k[ 02 ]- 
A similar approach is used to describe the local concentration of growth factors (e.g.. Vascular 
Endothelial Growth Factors) denoted by [VEGF]. The latter is responsible for the growth of new 
blood vessels, which is especially important for larger tumors. Initially, the tumor can benefit 
from the existing vasculature, but when the tumor occupies a larger volume the oxygen supply 
doesn’t satisfy anymore and the cells are obliged to use their ability to ask for new vessels, 
by secreting VEGF. Endothelial cells - the building blocks of blood vessels - react and move 
chemotactically towards the hypoxic regions. The corresponding reaction diffusion equation for 
VEGF reads: 


dt\VEGF]{x, t) = Dy'eg¥'^^\VEGF]{x, t) — ij)Y-EG¥nv{,x,t)\VEGF]{x,t) + 5'vEGF(a^, t) 

diffusion exchange with blood production 

— (5vEGFVEGF(a:,t) 

^ "V 

decay 

p Fit) 

>S'[VEGF] [x, t) = fc[vEGF] E E p(t)^( [^^^^int]i,p(0 [VEGFintjthr) 

p—1 i—1 


(14) 


Parameter 

Oxygen 

VEGF 

units 

D 

0.0014 

6 X 10-4 

cm^/min 


6 

6 X 10 "^ 

cm / min 

5 

0 

0.6 

min'^ 

ho2] 

-13 


min'^ 

fc[VEGF] 


0.6 

min'^ 

[02]ref 

20 


mmHg 


Table 2: Parameter values reaction diffusion equations 
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3 Variance reduction 


In this section, we propose a variance reduction algorithm similar to the technique used in 32 to 


simulate bacterial chemotaxis. The main differences are due to the fact that (i) the model is not 
conservative; and (ii) the internal dynamics only relates to cell division, apoptosis and VEGF 
secretion and not to advection-diffusive behavior. 

As in 32 , the algorithm relies on the combination of three simulations: a stochastic simulation 


with the full microscopic model, as well as with a coarse approximation, combined with a deter¬ 
ministic grid-based simulation of the coarse model. The full microscopic model uses an ensemble 
of Ip{t) particles with state variables: 


[p53],.p(t), [VEGFi„t]qp(i)}-^^i^ 


(15) 


As the coarse agent-based model, we conceptually consider an agent-based model in which the 
internal state has been suppressed and only the position remains: 


ipii) 




(16) 


So no internal dynamics is present, cells cannot divide, die or secrete VEGF. (In practice, we 
will use the results obtained from the full microscopic model, in which we neglect apoptosis and 
cell division, see later). The only dynamics is motion, which can be modeled with a PDE for the 
population density (see equation @). We call this coarse approximation the control process. 
We also introduce the formal semigroup notation: 




with Lpiup) = -DpV^ - XpV 


np(x,t) 1 - 


np{x,t) 


V[VEGF]{x,t) 


^p,max 


(17) 


that represents the exact solution of the macroscopic partial differential equation (101. In prac¬ 


tice, the solution will be approximated by a deterministic solution on a grid. It should be 
clear that the advection-diffusion behaviour in both agent-based models is identical. Thus, the 
only difference between the two models occurs when cells divide or die. Assuming no reac¬ 
tions take place, the three processes thus have the same expectation. This observation leads to 
the following variance reduction algorithm. As an initial condition, we start from /p(0) parti¬ 
cles sampled from specific probability densities, resulting in the number density np{x,0). For 
each particle, we choose a given internal state, for instance 'hi,p(0) = ^i,p(0) = 0, Ci,p(0 ) = Oj 
[ p53]i^p(0) = 0,[VEGFint]i.p(0) = 0, 1 < f < /p(0), I < p < P. (These internal states could 
also be sampled from an appropriate probability distribution.) Additionally, we introduce the 
variance reduced measure h{x,t), which we initialize as np{x,0) = np{x,0). We denote the time 
step St and the discrete time instances t^ = £St, £ = 0, 1, ... 

Algorithm 1 (Variance reduction for tumor growth). We advance the variance reduced number 
density n{x,t) from time t^ to t^'^^ as follows: 

• Evolve the particle states (151 from t^ to t^'^^ using the agent-based model §• 

• Compute the number density for the stochastic microscopic model using @ , as well as the 
number density for the coarse process as 




n(((a:,t^+^) = ^ Wi^p{t^)5x^^^(t^+^) 


(18) 


i=l 


i.e., we compute the number density for the control process based on particle positions 
and velocities at time t^^^, taking into account only the particles that were present in the 
simulation at time t^. 
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Evolve the control number density np{x,t) using a grid-based method based on ( |10[ ) and add 
the reactions (the difference in number density due to cell division and apoptosis) 


,{x,t^~^^) := np{x,t^) e^*^ np{x,t^^^) — np{x,t^~^^) 


(19) 


Next, we will prove that the proposed estimator for the population densities Up are unbiased 
and that the algorithm indeed reduces the variance on the mean population density. To this end, 
we define the so-called reaction field: 

Definition 2 (Reaction field). The control process differs from the full microscopic model in 
the way that there are no births, deaths or VEGF-secretion events. The direct influence on the 
population density can be summarized by the Reaction field defined as: 

Rp{x, = np{x, t*+^) — np{x, t*+^) 

/p(d+i) 7p(d) 

E{t‘) 


( 20 ) 


Definition 3 (Deterministic control density). ITe also introduce a shorthand notation for the 
control density calculated with the macroscopic evolution equation: 


n{x, ) := n{x, t)e‘ 


StL^ 


( 21 ) 


This procedure will be repeated after reinitializing the control density hp{x,t) = n{x,t). 
The importance of reinitialization can be illustrated by looking into the following hypothetical 
situation. Suppose the fth cell of type p divides at time t = t*, and hence cell Ip -I- 1 is born. 
At time t > t*, this newborn cell has moved randomly through the domain. Apart from this 
random motion, it also has influenced the environment along its track. Those events cannot be 
taken into account without reinitialization. 

Theorem 4 (Unbiased estimator). The algorithm described above yields an unbiased estimator 
for the population density Up. 

Proof. Assume that discretization errors are absent. Then, we can calculate the expectation 
value of rip based on equation (19) as follows: 

E [fipix, t’-'^^)] = E [np[x, -I- E [np{x, - E [np{x, 

By using the definition of sequentially the definition of hp (see equation ( [^ ) and the linearity 
of E, we can conclude that E[fip] = E[np] and hence hp is indeed an unbiased estimator. □ 


4 Results 

In this section, we will illustrate the performance of the variance algorithm described above 
with various numerical experiments. The cells are living on a 50 x 50 square grid. By default 
2000 normal cells are uniformly distributed over the whole domain. A small tumor consisting 
of 200 cancer cells are initially normally distributed with mean 0.25Aa: and standard deviation 
0.05Aa; in close to the left vessel. We simulate the system over 1920 timesteps (or 40 days). 
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The whole set of default parameters is summarized in table The normal tissue on the other 
hand is uniformly distributed over the whole domain. To initialize the agent-based simulation, 
we sampled /p(0) particles from the corresponding distribution. Remark that we mostly use a 
high number of particles to discretize the population density stochastically in order to make the 
agent-based model consistent with the continuum description. Hence, the equation to calculate 
the cell number density can be rewritten as: 

Ip{i) 

np{x,t) = ^ q^,pWi,p{t)5xi^^(t) (22) 

i=l 


where an additional weight qi^p is attached to each particle. This implies that each particle 
has a lower mass. The total mass is qi,pWi,pit). During the numerical experiments, we 

will use weights Qp independent of both the specific i-th particle of population p and of the time. 
Further, we initialize the environment as follows: two straight vessels at a; = 20Aa: and x = 40Ax, 
corresponding to a moderate vascular density of 50 cm^/cm^ (see 29 ). The latter results in av¬ 


erage oxygen concentrations, meaning that cells are proceeding through the cell cycle at a speed, 
which is slightly higher than half maximal. More details concerning realistic vascular densities 
and oxygen concentrations can be found in the supplementary material provided with 


29 


The macroscopic equations are simulated using a simple Euler discretization for the time deriva¬ 
tive and a second order central finite volume to discretize the spatial derivative. In the first three 
experiments, we have chosen for an explicit method. Further, the linear systems originating from 
the reaction-diffusion PDEs modeling the environment are solved using a conjugate gradient al¬ 
gorithm ( [^). The default choice of discretization, time-step and number of cells can be found 
in table 0 

The discussions corresponding to each of the individual experiments are organized as follows. 
First we consider the evolution of the population densities and the environment. Afterwards, 
we take a closer look at the variance with and without variance reduction. The results of the 
experiments are obtained by averaging out over 100 realizations. 


Remark 1 (Color code). During the numerical experiments, we adopt the following color code 
to describe the different quantities: 


• A colormap from white (low) towards gray (high) is used for the mean population densities 
(both normal and cancer). 

• A colormap from white (low) towards blue (high) is used for the variance on the mean 
cancer cell density (with and without variance reduction). 

• An additional colormap from green (negative) towards white(zero) and blue(high) is used 
to denote the covariance between the density and the reaction field. 

• A colormap from white (low) towards red (high) is used for the mean oxygen distribution. 

Remark 2 (units). We use minutes as default time unit and cm as the spatial unit. Those will 
be omitted in the figure titles for compactness. 


4.1 A small-scale experiment with a static vasculature 

Population densities In figure we have plotted the population density of the normal and 
cancer tissue, along with the oxygen concentration at time t = 5.76 x 10® min. 
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Parameter 

Normal 

Cancer 

EC 

units 

Dp 

0 

5 X 10-y 

1 X 10-3 

cm^/min 

4(0) 

2000 

200 

0 

#particles 

Qi,p 

1 

0.5 

1 

dimensionless 

St 

30 

30 

30 

min 

^p,max 

1 

2 

2 

#particles 

Ax 

4 X 10-3 

4 X 10-3 

4 X 10-3 

cm 

a 

0.5Aa; 

0.25Aa:: 


cm 

b 

0.5Aa; 

0.05Aa:: 


cm 


Table 3: Default parameter set used for the numerical experiments 



Q.l Q.l 0.1 


Figure 1: Mean cellular densities -normal tissue: (left panel), cancer cell density: (middle 
panel)- and mean oxygen concentration (right panel) calculated using variance reduction at 
time t = 5.76 x 10^ min. 


The tumor immediately influences the normal tissue in the sense that a significant amount of 
normal cells die in this cancerous region. They literally have to make place for the growing tumor. 
Besides cells also die due to lack off oxygen in between the two existing vessels. In the meantime, 
the tumor starts to grow along the vessel until the tissue is locally saturated, meaning that 
t) > ^p,max- The high birth rate can be explained by the high oxygen concentration. 
Furthermore, the cells also diffuse in the other directions due to the random Brownian motion. 
This evolution can also be seen as an illustration of the “go or grow”-paradigm, which is identified 
as an important characteristic of the aggressiveness of the tumor [O 
This process continues until the tissue is fully saturated along the leftmost vessel. Remark that 
none of the cells is able to cross the low oxygen zone. They would all die due to the hypoxic 
environment. 
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Evolution of the variance. We investigate how the population densities shown in hgure 
influence the variance. Figure illustrates the elements contributing to the variance (see equa¬ 
tion (251). We first consider the variance with and without reduction in more detail. Combining 
the definitions of both variance and rip, rip yields: 


VaT[np{x,t)] = VaT[np{x,t) + Rp{x,t)] (23) 

= Var[np(a;, t) + Yar[Rp{x, t)] -I- 2Cov{np{x, t), Rp{x, t)) (24) 

Var[fip] = Var[hp]-|-Var[i?p(x, t)]-I-2Cov(hp(a;, t), i?p) (25) 

the variance on the reaction field (left), the variance on the corresponding control densities - 
without variance reduction and with variance reduction - (middle) and the covariance 
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between the reactions and the control densities (right). We compare the results without (first 
row) and with variance reduction (second row). The variance on the reaction field is stretched 
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Figure 2: Evolution of the factors contributing - Var(i?2(a;, t)) (left), Var(n2(a;, t)) (middle) and 
Covar((i?2, ri 2 ){x, t)) (right) - to variance on mean cancer cell density with (second row per time) 
and without variance reduction at time t = 5.76 x 10® min. 

along the leftmost vessel, as is the tumor itself. To explain the pattern in more detail, we have 
to compare the variance on the reaction field with the corresponding density. The variance is 
especially high just next to the largest concentration of the tumor, where the concentration of 
reactions is high due to the combined effect of the relatively high number of cancer cells and the 
fact the tissue is not fully saturated yet. Remark that the variance on the reaction field does not 
depend on the variance reduction, since it is fully determined by the results of the agent-based 
simulation. 

However, the image is completely different for the variance on the corresponding control densities. 
Observing the middle picture on the first row leads to the conclusion that the noise is dominated 
by random motion since the variance is larger in the middle of the tumor than at the border where 
most of the reactions take place. In contrast, after applying variance reduction the variance is 
fully determined by the reactions, implying that the variance is mostly filtered by the algorithm. 
The above analysis (see figure]^ of the evolution of the variance clearly demonstrates the strong 
correlation between the oxygen concentration and the variance on the population densities. Again 
a more detailed view of the evolution can be found in the supporting material |6. 2 1 

To illustrate the performance of the algorithm in another way, we have taken some slices 
- at y = 0.04 cm, y = 0.1 cm and y = 0.14 cm respectively - of the cancer cell density at 
t = 5.76 X 10® min. In figure they are plotted along with their 95% confidence interval. The 
results based on the full stochastic model are plotted in red, while the results using the variance 
reduction algorithm are colored in blue. It can easily be seen that there is indeed a significant 
reduction and that the results using the variance reduction algorithm are consistent with the 
original results in the sense that they closely approximate the solution from the full stochastic 
model and that the variance is reduced significantly. 
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Figure 3: mean cancer cell density and reliability interval at time t = 5.76 x lO^min at y = 
0.04 cm (left), y = 0.1 cm (middle) and ?/ = 0.14 cm (right) with (blue) and without variance 
reduction (red). 

4.2 Experiment on a larger domain 

As pointed out before, our lattice-free approach allows to rescale the system in a straightforward 
way. Since the cost mainly depends on the number of particles and only marginal on the the 
domain size, it is possible to consider to perform a similar experiment on a rescaled (coarser) grid. 
To illustrate this we perform the simulation with Aa; = Ay = 1.26 x 10“^ cm, corresponding 
to a domain of 0.4cm^, corresponding with an upscaling of a factor 10. The normal tissue 
initially consists of 2 x 10"^ particles and a tumor of 1000 cells. In figure]^ we have plotted the 
evolution of the cellular distributions of the different cell types and the corresponding oxygen 
concentration. From the plot in the left column, one can see that normal cells are multiplying 
along the rightmost vessel since the left vessel is fully occupied by the tumor and there is not 
enough space for both the tumor and normal cells. In the rest of the domain the normal tissue 
is reduced to a minimal level due to lack of oxygen and the influence of the tumor. In figure 


(ni(x, t = 5.76 X 10®)> (n2(^, * = 5.76 X 10®)> (102](^, * = 5.76 X lO®)) 



Figure 4: Mean population densities and oxygen distribution in a large scale setting at time 
t = 5.76 X 10® min. 

we examine the influence of the variance reduction algorithm on the variance on the resulting 
tumor cell density as a function of time. Comparing the variance plot with (right panel) and 
without (left panel) give rise to the observation that the algorithm again yields a reduction of 
the variance both in the center and at the border of the tumor. This implies that the border of 
the tumor can be estimated in a more accurate, which determines the harshness of the tumor. 
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Figure 5: Evolution of the variance on the mean cancer cell density in a large-scale setting with 
and without variance reduction at time t = 5.76 x 10® min, zoomed in on the left vessel. 

4.3 Variance reduction for sprouting angiogenesis 

As a last experiment we will examine the performance of the algorithm in the case where the 
vasculature is also updated dynamically according to the model outlined in the section A 
small tumor mass of initially 100 cells - sampled from a uniform distribution, with parameters 
02 = O.SAx, 62 = O.lArc. The population density is discretized with 200 cells, i.e. qi ^2 = 0.5. 
Further, we have chosen ZI 2 = 1 x 10“® cm^/min. The other parameters are set to the default 
values outlined in table In contrast to the previous experiments, the cancer cells are now 
able to cause extension of the vascular network according to their needs. The resulting oxygen 
distribution reveals a strong correlation with the cancer cell distribution itself, meaning that 
the tumor is fully vascularized now and can grow further. A small fraction of the tumor even 
managed to reach the second vessel supported by some new branches in the vascular network 
created in response to the high VEGF gradients. Next, we investigate how the variance reduction 


t = 5.76 X 10®)> (n2(^, * = 5.76 X lO®)) (102](^, * = 5.76 X lO®)) 



Figure 6 : Mean cellular densities and oxygen distribution at time t = 5.76 x 10® min with 
dynamical vasculature 

algorithm is performing in this setting of dynamic vascularization. In figure we have plotted 
the variance on the mean cancer cell density with (n 2 {x,t)) variance reduction on the right and 
and without variance reduction on the left. 

A comparison of the two variance plots in figure leads to the conclusion that the variance 
is reduced everywhere. Along the leftmost vessel, the algorithm was even able to eliminate the 
noise completely. Indeed,the tissue is saturated here, so new cells are born here. Just outside 
this zone of maximal saturation the tissue is not so dense giving rise to more births and a higher 
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Figure 7: Variance with and without variance reduction on the mean cancer cell density with a 
dynamic vasculature at time t = 5.76 x 10^ min 


level of noise here. In this region, the noise is proportional to the local density itself. 


4.4 Fast diffusing cancer cells 


Motivated by the hypothesis that the diffusion coefficient can be related to the aggressiveness 
of the tumor, we investigate the situation where cancer cells have a higher diffusion coeffi¬ 
cient. Swanson et al. have shown 35 that glioma’s with a higher diffusion coefficient have 
a higher probability to cause metastases, which is obviously an important characteristic for 
the long-term survival probability of the patient. Apart from the modified diffusion coefficient 
Z ?2 = 5 X 10“^ cm^/min, we adopt the same initial configuration as in the previous experiment. 
Remark that the simulations are performed with a smaller time-step {5t = 0.3 min) in order to 
fulfill the CFL-condition corresponding to the macroscopic equation. 


Populations In figure]^ the evolution of both the normal tissue and the tumor are shown along 
with the local oxygen concentration at t = 1.152 x 10^ min. As in the previous experiment, the 
normal tissue density is the result of the cell deaths due to the presence of the tumor, while 
the normal cells are more sensitive to hypoxic environment. The tumor, on the other hand, has 
diffused through the normal tissue significantly on this short timescale, without consuming too 
much oxygen. Indeed, the cancer cells have already covered a large distance within a rather 
short time interval, meaning that the tumor exhibits the 50 -phenotype, rather than the grow- 
phenotype as it was the case in the first experiment. Despite the fact that the tumor didn’t 
cause a lot of damage, it is potentially dangerous since it can stay more or less invisible for a 
long time and as soon as the tumor reaches a vessel it is possible that cells invade a vessel and 
give rise to metastatic spread of the cancer. 


Evolution of variance As before, we also examine the variance on the mean cancer cell 
density with and without variance reduction. Without variance reduction, the resulting variance 
is proportional to the density itself, suggesting that the variance is mainly caused by the random 
jumps. Obviously, the latter will be higher in zones with more cells. When variance reduction 
is applied, the variance is reduced with at least a factor 100 pointwise and moreover the plot 
reveals a clear pattern, which is again related to the oxygen concentration. 
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Figure 8: Evolution of the cellular populations as a function of time (at time t = 1.52 x 10"^ min, 
with fast diffusing cancer cells. 
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Figure 9: Evolution of the variance on the mean cancer cell density with and without variance 
reduction at time t = 1152 min. 

5 Discussion 


We developed a novel variance reduction technique specihcally suited to reduce the noise of 
agent-based models with birth and death events, as it is the case in our model for tumor growth. 
We proved that the algorithm outlined in section gave rise to an unbiased estimator and the 
variance is determined by the births and deaths. The performance was illustrated numerically 
in different possible regimes characterizing different aspects of tumor growth such as sprouting 
angiogenesis, highly diffusive cancer cells and large-scale systems. The proposed algorithm is 
based on the idea of control variates, since the evolution of the system without reactions is 


known deterministically via the macroscopic equation (10). 


A valuable extension would be to combine this algorithm with other variance reduction techniques 
such as importance sampling. It is self-evident that an accurate and efficient simulation of all 
the different aspects of the system is crucially important for the reliability of the system as 
a whole. Apart from that, we will also extend our model with important features such as 
haptotaxis in response to the extra-cellular matrix and include a more sophisticated model for 
stemcellness [6l[7l[26l since it was identified as one of the hallmarks of cancer [TtI . Another track 


worthwhile further investigation is to apply our technique to patient-specific data. For instance, 
patient-specific data, like MRI-images or blood parameters, could be used as a specific initial 
configuration 25 . Finally, this algorithm can also be applied on related systems such as bone 


fracture healing and other application where we are interested in macroscopic behavior, but with 
agent-based features characterizing the dynamics. 
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6 Supporting Information 

6.1 SI Video 

Evolution of the population densities in the small scale setting. Evolution of the 
mean normal and cancer cell density from time t = 0 till time t = 5.76 x 10® min. The initial 
configuration corresponds with the first numerical experiment. 

6.2 S2 Video 

Evolution of the variance in the small scale setting. Evolution of the variance on the 
mean normal and cancer cell density from time t = 0 till time t = 5.76 x 10® min. 
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